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Abstract 

Most mathematical models of collective cell spreading make the standard assumption that the cell diffusivity and cell 
proliferation rate are constants that do not vary across the cell population. Here we present a combined experimental and 
mathematical modeling study which aims to investigate how differences in the cell diffusivity and cell proliferation rate 
amongst a population of cells can impact the collective behavior of the population. We present data from a three- 
dimensional transwell migration assay that suggests that the cell diffusivity of some groups of cells within the population 
can be as much as three times higher than the cell diffusivity of other groups of cells within the population. Using this 
information, we explore the consequences of explicitly representing this variability in a mathematical model of a scratch 
assay where we treat the total population of cells as two, possibly distinct, subpopulations. Our results show that when we 
make the standard assumption that all cells within the population behave identically we observe the formation of moving 
fronts of cells where both subpopulations are well-mixed and indistinguishable. In contrast, when we consider the same 
system where the two subpopulations are distinct, we observe a very different outcome where the spreading population 
becomes spatially organized with the more motile subpopulation dominating at the leading edge while the less motile 
subpopulation is practically absent from the leading edge. These modeling predictions are consistent with previous 
experimental observations and suggest that standard mathematical approaches, where we treat the cell diffusivity and cell 
proliferation rate as constants, might not be appropriate. 
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Introduction 

Collective cell spreading plays an important role in development 
[1], repair [2-5] and disease [6]. One way of improving our 
understanding of the mechanisms that influence collective cell 
spreading is to develop and implement a mathematical model that 
can both mimic existing experimental observations as well as 
suggesting new experimental options for studying collective cell 
spreading [7]. Such mathematical models have provided key 
insights into several biological systems. For example, Greenspan's 
model [8] of tumor growth provided a potential explanation of the 
observed spatial structure in tumor spheroids, while Gatenby and 
Gawlinski's model of tumor spreading into surrounding tissue [9] 
predicted the formation of a gap between the two types of tissue 
that was later verified experimentally [7] . 

Almost all mathematical models of collective cell spreading 
processes make the simplifying assumption that the population of 
cells can be treated as a uniform population. For example, Maini 
and coworkers [2,3] studied a scratch assay and showed that the 
solution of a reaction-diffusion partial differential equation led to 
constant-speed, constant-shape moving fronts that were consistent 
with experimental measurements. Similarly, Sengers and cowork- 
ers [10, 1 1] studied a circular cell spreading assay and showed that 
the solutions of an axisymmetric reaction-diffusion equation 
matched the time evolution of the observed experimental cell 
density profiles. These studies made an implicit assumption that 
the motion of cells within the population could be described using 
a constant value of the cell diffusivity D, and that the proliferation 



rate of cells could be described by a constant value of the cell 
proliferation rate, X. Similar assumptions are often made in 
discrete models of collective cell motion [12]. For example, Cai 
and coworkers [13] used a random walk model to study 
experimental observations of a scratch assay where the motility 
of isolated individual agents and the birth rate of isolated 
individual agents in the discrete models were treated as constants. 
Similarly, Binder and coworkers [14] applied a discrete random 
walk model of cell migration and cell proliferation on a growing 
tissue while Khain and coworkers [15] applied a discrete random 
walk model incorporating cell migration, cell proliferation and 
cell-to-cell adhesion to a scratch assay performed with glioma cells. 
Khain's discrete model treated the cell motility, cell proliferation 
rate and cell-to-cell adhesion strength as a constant for each 
isolated agent in the simulations. 

In contrast to many mathematical models, there are a range of 
experimental observations which suggest that cell motility and cell 
proliferation rates are not constant and might vary considerably 
amongst a population of cells. For example, during the develop- 
ment of the drosophila nervous system, time-lapse observation of 
individual glia cell migration and proliferation have reported the 
formation of glial chains which appear to be an essential 
component of normal development [16,17]. Time-lapse imaging 
and cell ablation experiments suggest that a certain subpopulation 
of the glial cells act as pioneer (or leader) cells, and that these pioneer 
cells guide the behavior of the remaining follower cells. A similar 
chain migration model has been proposed to explain time-lapse 
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observations of the development of the enteric nervous system 
which involves a population of precursor cells, called neural crest 
cells, moving along the developing intestines in the form of chains 
of cells [18-22]. The details of this developmental system have 
been studied experimentally and the results suggest that cells at the 
leading edge of the population follow directed trajectories whereas 
cells located behind the leading edge of the population followed 
less directed, more random trajectories [19]. These observations 
have been recently incorporated into a discrete mathematical 
model of observed behavior in a related experimental system 
[23,24] where it was found necessary to make an explicit 
distinction between the behavior of pioneer and follower cells to 
replicate the observed patterns [25]. 

Experimental observations that are consistent with the existence 
of pioneer and follower cells have also been made in various in vitro 
assays. For example, Cai and coworkers recorded trajectories of 
individual cells within a scratch assay and showed that cells at the 
leading edge of the population moved along trajectories that were 
qualitatively different to other cells located behind the leading 
edge [13]. Distinct roles for pioneer and follower cells have been 
observed in cell populations that interact with collagen fibres [26] 
and in two-dimensional monolayers of cells that have been 
wounded [27]. Other biological systems which suggest a role for 
pioneer and follower subpopulations of cells include the immune 
system [28-30], three — dimensional tumor spheroid growth [31] 
and various aspects of development [32,33]. We note that, very 
recently, heterogeneity amongst circulating tumor cells in patients 
with advanced primary cancer has been proposed to explain 
variations in metastatic disease patterns [34]. 

In this work we investigate whether an apparendy homogeneous 
population of motile cells is composed of functionally distinct 
subpopulations that could be interpreted as a pioneer subpopu- 
lation and a follower subpopulation. This investigation makes use 
of both experimental measurements as well as a simplified 
mathematical model of collective cell behaviour that we use to 
represent both individual cell behavior and the emergent collective 
behavior of the entire cell population. We perform a three- 
dimensional transwell assay [35] where we stop the experiment 
after a relatively short period of time and remove those cells which 
have moved through the porous membrane as well as those cells 
which have not moved through the porous membrane. Both these 
populations of cells are cultured separately, and individual cell 
trajectories are recorded so that we can investigate whether there 
are any differences between the two groups of cells. Our 
experimental measurements are interpreted using a discrete 
three-dimensional mathematical model of cell migration in a 
transwell. Although, in principle, our mathematical model can be 
used to study a very general population of cells where each cell has 
a unique motility and proliferation rate, we take the simplest 
possible approach and interpret our experiments by making the 
assumption that the total population is composed of just two 
subpopulations which we refer to as (i) subpopulation 1 with cell 
diffusivity, D\, and cell proliferation rate, %\, and (ii) subpopula- 
tion 2 with cell diffusivity, Z>2, and cell proliferation rate, X2. Using 
our model we show that our transwell results are consistent with 
the hypothesis that the two subpopulations are distinct since we 
find D\ >Z>2. Although we make no experimental measurements 
of collective behavior involving cell proliferation, we conclude by 
presenting some simulations of a scratch assay where proliferation 
plays an important role. In these simulations we treat the entire 
population as two interacting subpopulations and our modeling 
suggests that an initially well-mixed population of cells can form a 
spatially organized spreading front of cells where the more motile 
subpopulation dominate at the leading edge of the spreading 



population whereas the less motile subpopulation is practically 
absent from the leading edge. 

Materials and Methods 

1.1 Experimental methods 

Mouse fibroblast feeder cells [36] (3T3 cells) (ATCC, CCL-92, 
Manassas, VA, USA) were used to perform the transwell migration 
assay. The 3T3 cells were cultured in Dulbeccos modified Eagle 
medium (DMEM; Invitrogen, Australia) supplemented with 5% 
foetal calf serum (FCS; Hyclone, New Zealand), 2 mM L- 
glutamine (Invitrogen) and 1% v/v Penicillin/ Streptomycin 
(Invitrogen) in 5% C0 2 and 95% air at 37°C. 

A schematic of the transwell apparatus is shown in Fig. 1A, and 
the assay was performed as previously described [37]. In brief, the 
3T3 cells were serum starved for four hours by incubating in 
serum free medium (SFM). The SFM was DMEM without FCS. 
The cells were harvested, and the flasks washed with phosphate- 
buffered saline (PBS; Invitrogen) followed by exposure to 0.05% 
trypsin-EDTA (Invitrogen) for one-to-two minutes at room 
temperature. The cell suspension was collected in a 50 mL falcon 
tube and centrifuged twice at 1000 rpm for five minutes to 
eliminate the trypsin. The supernatant was discarded and the 
pellet re-suspended in 10 mL of SFM. The viable cells were 
counted using a trypan blue exclusion test and a haemocytometer. 

Fifty thousand 3T3 cells suspended in SFM were seeded into the 
upper compartment of a 12 |xm pore transwell (Corning, New 
York, USA) where the under-surface of the porous membrane had 
been pre-coated with with 10% FCS. Each transwell was placed in 
a 12 — well plate which was incubated for two hours at 37°C with 
5% C0 2 and 95% air. After two hours, those cells that had moved 
into the lower compartment and those cells that remained in the 
upper compartment were collected separately using 0.05% 
trypsin-EDTA. The transwell inserts were first rinsed with PBS 
and then trypsin was introduced into the upper and lower 
compartments to collect the two groups of cells separately. The 
collected cells were centrifuged at 1000 rpm for five minutes to 
remove trypsin and re-suspended in 1 mL of 3T3 medium. 

Both groups of cells were separately re-seeded onto a 24-well 
tissue culture plate and monitored using a widefield microscope 
(Leica, Australia). Images were captured at five minute intervals 
over a period of 16 hours. 

1.2 Mathematical modeling tools 

An interacting random walk model, that explicitly incorporates 
cell-to-cell crowding effects, is used to simulate the experiments. The 
model is implemented on a three — dimensional square lattice with 
spacing A. Each site is indexed (ij,k), where i, j, kel, and has 
position (x,y,z) = (iAjA,kA). A random sequential update method 
[38] is used to perform the simulations so that if there are N(t) 
agents at time t, during the next time step of duration 1, N( t) agents 
are selected at random, one at a time, and given the opportunity 
move with probability P^e[0,l], where /= 1,2,3, . . . ,N(t). Speci- 
fying different values of P l m allows different agents in the model to 
move with a different, unique, motility rate. A motile agent at site 
(ij,k) attempts to step to: (i) (i± lj,k) with probability (1 ±p x )/6, 
(ii) (ij+\,k) with probability (l+p y )/6, or (iii) (ij,k+l) with 
probability (1 ±p z )/6. The parameters p x e[ — 1,1], P y e { — 1>1] and 
p z e[ —1,1] control the motility bias: setting p x = p v = p z = 0 means 
that the motion is unbiased. If an agent attempts to step to an 
occupied site, then that motility event is aborted. Once the N(i) 
potential motility events have been assessed, another N(t) agents are 
selected at random, one at a time, and given the opportunity to 



PLOS ONE I www.plosone.org 



2 



January 2014 | Volume 9 | Issue 1 | e85488 



Do Pioneer Cells Exist? 



Upper compartment 

Porous membrane 
Lower compartment 




D 

15 



20000 

<Ni> 

<N 2 > 10000 



200 400 
D (|um 2 /minute) 



D (|um /minute) 




20000 



10000 



100t 200t 
time (t) 



300t 




100t 200t 
time (x) 



300x 



Figure 1. Experimental results and three dimensional mathematical modeling results for a transwell assay. Schematic of a transwell 
apparatus illustrating that the cylindrical insert is 12 mm in diameter, and that the upper and lower compartments are separated by a porous 
membrane (A). At the conclusion of the two hour transwell migration assay those cells that moved into the lower compartment were collected and 
placed on a cell culture plate. The trajectories of individual cells were recorded over a period of 1 6 hours. The white scale bar is 1 00 |im (B). Similarly, 
at the conclusion of the two hour transwell assay those cells that remained in the upper compartment of were collected and placed on a tissue 
culture plate. The trajectories of individual cells were recorded over a period of 16 hours. The white scale bar is 100 |im (C). The trajectories of 20 
individual cells from those that moved into the lower compartment were analyzed to produce 40 estimates of the cell diffusivity D, shown as a 
histogram (D). The average cell diffusivity of those cells that had moved into the lower compartment was <£>> = 102 |im 2 /minute. The trajectories of 
20 individual cells from those cells that remained in the upper compartment of the transwell were analyzed to produce 40 estimates of the cell 
diffusivity D, shown as a histogram (E). The average cell diffusivity of those cells that had not moved into the lower compartment of the transwell was 
<£>> = 31 |im 2 /minute. Three dimensional simulation results of a transwell assay initialized with 25000 cells from subpopulation 1 and 25000 cells 
from subpopulation 2 (F-G). Simulation results show <iVi> and <W2>, corresponding to the average number of cells associated with subpopulation 1 
and subpopulation 2 remaining in the upper compartment as a function of time. The average simulation results were obtained using M=100 
identically prepared realizations of the three dimensional random walk model. Simulation results correspond to two cases: (i) identical 
subpopulations with P) n = P] n = 1 and Pj = Pj=0 (F), and (ii) distinct subpopulations with P)„ = \, ^=0.3 and Pj = .Pj| = 0. 
doi:1 0.1 371 /journal. pone.0085488.g001 



proliferate with probability Pj!,e[0,l]. We model proliferation with 
an unbiased mechanism whereby a proliferative agent at (ij,k) 
attempts to deposit a daughter agent at (/ + 1 J,k), (ij + 1 ,k) or 
(ij,k+ 1), with each target site chosen with equal probability 1/6. 
Potential proliferation events that would place an agent on an 
occupied site are aborted [35,39,40]. 



This basic modeling framework will be applied to two different 
experimental scenarios. First, we will apply this three-dimensional 
model directly to the geometry of the transwell apparatus as we 
have done previously [35]. Second, we will consider a simpler 
two-dimensional application of the model which is consistent with 
a two-dimensional in vitro assay, such as a scratch assay [41]. We 
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note that the three-dimensional modeling framework can be used 
to simulate a two-dimensional assay simply by considering a 
three-dimensional lattice with a single layer in the vertical 
direction, so that k = 1 . In the simpler two-dimensional format a 
motile agent at site (ij) will attempt to step to: (i) (i+lj) with 
probability (1 + p x )/4, or (ii) (ij + 1) with probability (1 ±p y )/4. 
Similarly, a proliferative agent at (ij) attempts to deposit a 
daughter agent at (;'+ Ij) or (ij + 1), with each target site chosen 
with equal probability 1/4. 

Although, in principle, our discrete modeling framework can be 
applied to a very general system by allowing every single agent 
within the population to have a unique motility and proliferation 
rate, we will implement our model using the simplest possible way 
to investigate the role of variability within the total population by 
making the assumption that the population is composed of two 
subpopulations: (i) subpopulation 1, which is composed of cells 
which have a probability of motility per time step of P x m and a 
probability of proliferation per time step of P p , and (ii) 
subpopulation 2, which is composed of cells which have a 
probability of motility per time step of P 2 m and a probability of 
proliferation per time step of P 2 p [40] . 

We would like to point out that while our mathematical model 
explicitiy incorporates physical interactions between cells in the 
population by incorporating cell-to-cell crowding and volume 
exclusion effects, our mathematical model is an idealization of 
collective cell behaviour. One important aspect that our model 
neglects is any consideration of biochemical interactions amongst 
the population of cells, which can play an important role in 
collective cell behaviour [20,21]. The neglect of such biochemical 
interactions is a standard assumption made in many mathematical 
modelling studies of collective cell migration [2-5,8,10,11,13,15] 
and the focus of our present work is not to build a mathematical 
model which incorporates every detail of collective cell migration. 
Instead, the focus of our present work is to investigate the role of 
variability amongst a population of cells since traditional mathe- 
matical models of collective cell behaviour routinely treat the 
motility of cells as a simple constant value across a population of 
cells [2,3,10,11,13]. Similarly, most traditional mathematical 
models of collective cell behaviour routinely treat the proliferation 
rate of cells as a simple constant value across a population of cells 
[2,3, 1 0, 1 1 , 1 3] . The aim of our work is to explore the validity of such 
assumptions and to use a simplified mathematical model to 
demonstrate the implications of such assumptions. 

Results 

2.1 Transwell Results 

2.1.1 Estimating the cell diffusivity. Once the cells were 
harvested at the conclusion of the two hour migration period in the 
transwell apparatus, those cells that had migrated into the lower 
compartment of the transwell (Fig. 1A) were collected separately 
from those cells that remained in the upper compartment (Fig. 1 A). 
These two groups of cells were placed on separate culture plates 
and individual cells within the two groups were imaged using time- 
lapse microscopy for a period of 16 hours so that we could 
characterize the motility of both populations. At the conclusion of 
the 16 hour period the time-lapse images were analyzed using 
ImageJ to record the trajectories of individual cells within the 
population [42]. For simplicity we will refer to those cells that 
migrated through into the lower compartment of the transwell as 
subpopulation 1, and those cells that remained in the upper 
compartment of the transwell as subpopulation 2. 



To characterize the motility we estimate the squared displace- 
ment for the x-coordinate and j-coordinate of each trajectory 

x 2 (t) = (x(t)-x(Q)f, y 2 (t) = (y(t)-y(0)) 2 , (1) 

where (x(t),y(t)) are the two-dimensional Cartesian coordinates 
of the cell after time t. An estimate of the random motility 
coefficient, also known as the cell diffusivity, in each orthogonal 
direction is then obtained by fitting a least-squares straight line to 
the data [43], 

x 2 (t) = 2D x t, y 2 (t) = 2iyt, (2) 

where D x and D y are the diffusivities in the x and y directions. We 
analyzed 20 randomly chosen cell trajectories from each subpopula- 
tion, being careful that we only considered trajectories that did not 
collide with other cells during the 16 hour observation period. This 
gave us 20 estimates of D x and D y for both subpopulations. Averaging 
these data, for both subpopulations, indicated that <D X >»<Z) J '> 
which is reasonable since the substrate is isotropic [44]. Therefore, for 
each subpopulation we pooled the D x and D y data which are 
presented as histograms in Fig. ID and Fig. IE for subpopulations 1 
and 2, respectively. In both histograms the data shows that the 
majority of observed trajectories are associated with a diffusivity in the 
range D< 100 urn /minute. However, both subpopulations con- 
tained some cells that were much more motile, and we observed 
some trajectories corresponding to cell diffusivity estimates with 
D x 500 urn 2 / minute. Averaging the 40 diffusivity estimates for each 
subpopulation gives <Z>i>«102 umVminute and <Z> 2 >«31 um 2 / 
minute. These results indicate that subpopulation 1 is, on average, 
approximately 3.3 times more motile than subpopulation 2. 

2.1.2 Discrete simulations using the transwell data. To 
investigate how the variations within the cell population could 
affect our interpretation of a transwell assay we apply the three- 
dimensional mathematical model to the transwell apparatus using 
the same procedure outlined previously in [35]. In brief, the 
transwell is cylindrical with an inner diameter of 12 mm and the 
3T3 cells are, on average, approximately 25 |J,m in diameter [45]. 
We represent the upper compartment using a three dimensional 
lattice with A = 25 urn. The three dimensional lattice has five 
layers in the vertical direction giving 1 <k< 5, and each layer is a 
square with length l</<480 and width l</<480. The length 
and width are chosen to accommodate the 25 urn 3T3 cells in the 
12 mm diameter transwell so that we have 12/(0.025) = 480. To 
represent the cylindrical geometry, all sites in the region 
(z-240) 2 + (/'-240) 2 <240 2 are labeled as active sites, meaning 
that they can can be occupied by agents. The remaining sites 
where (/ - 240) 2 + (/ - 240) 2 > 240 2 are labeled inactive sites, which 
cannot be occupied by agents. Each layer in the lattice contains 
7i240 2 x 180956 active sites so that our model can accommodate 
up to 5(7i240 2 ) x 904778 agents. The porous membrane separates 
the upper and lower compartments and is approximately 15% 
pore space [35]. To model the porous membrane we randomly 
select 15% of the active sites on the lower (k= 1) layer of the lattice 
and assume that these sites, called downward permeable sites, represent 
a pore in the membrane. The remaining 85% of active sites on the 
lower (k = 1 ) layer are downward impermeable sites. In our model a 
motile agent residing on a downward impermeable site (ij,k), steps to 
(i) (i+lj,k) with probability (l+p x )/6, (ii) (ij+\,k) with 
probability (\+p y )/6, and (iii) (ij,k+l) with probability 
(l+p z )/6 and (ij,k—\) with probability zero owing to the 
presence of the porous membrane. In comparison, a motile agent 
residing on a downward permeable site (ij,k) is permitted to move in 
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the negative z direction in the usual way as this agent is not 
blocked by the membrane. 

During a transwell assay cells are placed in the upper compartment 
and rapidly settle onto the porous membrane [35] . We model this by 
placing agents on the lattice to mimic the way that cells are distributed 
after they have setded onto the membrane. For example, to model 
our experiments described in Section 1.1 we initially randomly 
occupy 27.7%= 100x50000/180956 of active lattice sites on the 
lower (k= 1) layer of the lattice. To represent the movement of cells 
in the transwell experiments we set p z = — 1 to prevent agents 
moving vertically upward which is consistent with our observations of 
cell movement in a transwell [35]. We also set p x = p y = 0, which is 
appropriate because we do not expect any bias in the horizonal plane. 
During the simulations some agents move vertically down through 
the pore space and we assume that these agents leave the system and 
no longer interact with other agents during that simulation. Any 
potential motility event that would place an agent on an inactive site, 
or on a site that is already occupied, is aborted [35,39,40] . Our model 
predictions are made by counting the number of agents leaving the 
system through the lower layer of the lattice. Since the algorithm is 
stochastic we present results by averaging over many identically 
prepared realizations of each simulation. 

Results in Fig. 1 F correspond to a simulation where the transwell 
experiment was initialized with 25000 agents from subpopulation 1 
and 25000 agents from subpopulation 2. In this case we make the 
standard assumption that both subpopulations are identical with 
= P 2 ^ = 1 . We note that many transwell assays are performed for 
periods of time that are much shorter than the cell cycle time [35]. 
This means that any increase in cell number due to cell proliferation 
is negligible during such experiments. To make our modeling 
consistent with this we set P p =P 2 = Q. Averaged modeling results in 
Fig. IF show the number of agents in each subpopulation that 
remain in the upper compartment as a function of time and we see 
that the time taken for both subpopulations to exit the upper 
compartment are the same. After approximately 100 time steps 
almost all of the agents have moved into the lower compartment. 
This result makes sense intuitively since we have specified that both 
subpopulations behave identically so we might have anticipated that 
both subpopulations will exit the upper compartment of the transwell 
at the same rate. We would like to point out that the results in Fig. IF 
are reported for an arbitrary duration of each time step, T. If, for 
example, we chose T=1.02 minutes, our simulations would 
correspond to D\ = D 2 « 102 u.m 2 /minute since we have 
A = (P> 2 )/(6t) and D 2 = (P 2 m A 2 )/(6t). 

Results in Fig. 1G correspond to a simulation where the 
transwell experiment was initialized with 25000 agents from 
subpopulation 1 and 25000 agents from subpopulation 2. In this 
case we assume that the subpopulations are distinct and we choose 
the motility parameters to reflect the differences we observed in 
the experimental data reported in Section 2.1.1. By choosing 
P x „=\ and P 2 m =0.3, we simulate two distinct subpopulations 
where subpopulation 1 is approximately 3.3 times more motile 
than subpopulation 2. Again, to be consistent with standard 
transwell protocols, we neglect any increase in cell number by cell 
proliferation by setting P l = P 2 p = 0 [35] . The averaged modeling 
results in Fig. 1G show that we observe very different behavior 
from the results in Fig. IF where we made the standard 
assumption that all the cells agents in the system behaved 
identically. In this case our modeling shows that subpopulation 1 
moves into the lower compartment much faster than subpopula- 
tion 2 (Fig. 1G). In particular, we see that after 100 time steps 
almost all of subpopulation 1 has moved into the lower 
compartment whereas almost 300 time steps are required for 



almost all of subpopulation 2 to move into the lower compartment 
(Fig. 1G). This difference in the behavior of the two subpopula- 
tions is expected since we have P x m > P 2 m , and so we anticipate that 
agents from subpopulation 1 are able to migrate around in the 
transwell much more efficiently than members of subpopulation 2. 
This would mean that agents belonging to subpopulation 1 are 
more likely to find the location of the pores in the membrane 
through which they can move into the lower compartment. We 
also note that the results in Fig. 1G are reported for an arbitrary 
duration of each time step t. If, for example, we chose t = 1 .02 
minutes, then this would correspond to D\ x 102 urn / minute and 
Z>2~31 pjnVminute which is consistent with our cell diffusivity 
estimates from our experiments as reported in Fig. ID and Fig. IE. 

In summary, our modeling results indicate that our interpretation 
of transwell assays could be very sensitive to differences amongst the 
motility rates of the cells. Examining the results in Fig. 1G indicates 
that if we stopped the simulation after a relatively short period of 
time, say 50t, then almost all of subpopulation 1 would have moved 
into the lower compartment while the majority of subpopulation 2 
would remain in the upper compartment. These averaged 
simulation results are consistent with our experimental observations 
in Fig. ID and Fig. IE since our experimental data indicates that the 
group of cells that moved into the lower compartment after a 
relatively short time period were, on average, more motile than the 
group of cells remaining in the upper compartment. 

2.2 Scratch assay 

Since our modeling results in Fig. IF and Fig. 1G imply that a 
transwell assay could be very sensitive to differences amongst the 
motility rate of the cell population, we now extend these ideas to a 
scratch assay [2,3,41]. Scratch assays are often performed in a 
narrow channel geometry where a confluent population of cells is 
wounded, or scratched, to reveal a sharp front that separates the 
confluent region from a vacant region. Typically, a scratch assay is 
monitored by measuring the location of the leading edge of the 
population as it spreads and the initially vacant region becomes 
occupied [2,3,41]. To model this we apply the discrete mathe- 
matical model on a two-dimensional lattice where each site is 
indexed (ij), and each site has position (x,y) = (iAJA). Here we 
choose A = 25 Mm to correspond to the diameter of 3T3 cells. We 
apply this model on a two-dimensional domain with 0<x<25 
mm and 0 < y < 1 .25 mm, to mimic the narrow channel geometry. 
Reflecting boundary conditions are applied along all boundaries. 

To be consistent with our results in Fig. IF and Fig. 1G, we 
consider the initial population of agents to be composed of two 
subpopulations. Each simulation is initialized so that the central 
region of the lattice, where 12<x< 13 mm, contains a confluent 
monolayer. This initial confluent monolayer contains, on average, 
50% of agents from subpopulation 1 and 50% of agents from 
subpopulation 2. Two different types of simulations are performed. 
In the first simulation (Fig. 2A-D) we make the standard 
assumption that both subpopulations are identical with 
Pl=P 2 m = l and P p = P 2 = 0.001. Unlike transwell assays, many 
scratch assays are reported for a period of time that is longer than 
the cell cycle time so that proliferation plays an important role 
[2,3,41] and therefore we include proliferation in these simulations 
[39]. Results in Fig. 2A-D show snapshots of the simulation after 
0, 1000, 5000 and 10000 time steps, where each time step has a 
duration T. These simulations show that the population spreads 
into the initially vacant region. Individual agent motility and 
proliferation events lead to the formation of two fronts, one 
moving in the positive x-direction and the other moving in the 
negative x-direction. The formation of such fronts is consistent 
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Figure 2. Two-dimensional modeling results for a scratch assay. Discrete snapshots of a two dimensional scratch assay in a narrow channel 
geometry with 0<x<25 mm and 0<y<1.25 mm (A-D, l-L). The initial condition for two different simulations corresponds to a confluent 
monolayer of agents in the central region of the domain, where \2<x< 13 mm. The initial population is made up of 50% subpopulation 1 (red disks) 



and 50% subpopulation 2 (blue disks). The first simulation corresponds to identical subpopulations with P]„ = P^„ = 1 and P p =Pi = 0.001 (A-D) and 

Pl=Pj =0.001 (l-L). Snapshots are shown after 0 (A,l), 1000 



the second simulation corresponds to distinct subpopulations with P'„ = 1, Pf„=03, 
(B,J), 5000 (C,K) and 10000 (D,L) time steps, where each time step has a duration of t. Both types of discrete simulation were repeated using M = 100 
identically prepared realizations to to produce the averaged density profiles for the case where both subpopulation are identical (E-H) and where the 
subpopulations are distinct (M-P). The numerical solution of Equation (5) was obtained for the initial condition given by Equations (8)-(9) and 
superimposed on the averaged discrete results (E-H, M-P). The numerical solutions were obtained using <5.y = 2.5 x 10~ 3 mm and <5? = 0.1t. 
doi:1 0.1 371 /journal.pone.0085488.g002 



with experimental observations where these fronts often move with 
constant speed [2,3] . We observe that the total population grows 
rapidly with time, and our simulation indicates that the two 
subpopulations remain well-mixed for all time and at all locations. 
The results in Fig. 2A-D are reported for an arbitrary duration of 
each time step, x. If, for example, we chose x = 1.53 minutes, this 
would correspond to D[=D2»W2 umVminute since we have 
D x =(P, 1 „A 2 )/(4t) and D 2 = (P 2 ,A 2 )/(4t) in two-dimensions. 
Similarly, choosing T=1.53 minutes corresponds to Xi=%2~ 
0.00065/minute since we have X\=P l p /x and X 2 = PjJx. This 
proliferation rate corresponds to a doubling time of tj « 18 hours 
since we have / f / = log e 2/X [35]. 

In the second simulation (Fig. 2I-L) we allow the two subpopu- 
lations to behave differently by setting P m = 1 and P 2 , = 0.3, so that 
subpopulation 1 is approximately 3.3 times more motile than 
subpopulation 2. Again, this difference in the motility rate between 
the two subpopulations is consistent with our experimental results in 
Fig. ID and Fig. IE. Since we have not made any measurements of 
the proliferation rate of cells we assume that both subpopulations 
proliferate at the same rate with Pj, = P 2 = 0.00 1 . Results in Fig. 21— 



L show a snapshot of the simulation after 0, 1000, 5000 and 10000 
time steps, where each time step has a duration of x. Just like the 
uniform population in Fig. 2A-D, we see that the population spreads 
into the initially vacant region of the domain and the model predicts 
the formation of two fronts, one moving in the positive A"-direction 
and the other moving in the negative x-direction. Again, the total 
population grows rapidly with time, however in this case our results 
indicate that the two initially well-mixed subpopulations remain 
mixed for a short period of time only (Fig. 2I-J) before becoming 
segregated at later times (Fig. 2K-L) where we see that the leading 
edge of the population is dominated by subpopulation 1 . This result 
implies that the leading edge of the spreading population becomes 
dominated by the subpopulation that is more motile. The results in 
Fig. 2I-J are reported for an arbitrary duration of each time step, x. 
If, for example, we chose x = 1 .53 minutes, this would correspond to 
D\ x 1 02 |im / minute and D 2 « 3 1 umV minute which is consistent 
with our experimental observations in Fig. 1D-E. 

2.3 Continuum description 

The simulation results in Fig. 2A-D and Fig. 2I-L correspond 
to single realizations of the discrete model. To provide additional 
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information about these simulations we consider M identically 
prepared realizations and generate averaged density profiles. In 
the m th identically prepared realization of the model, site (ij) can 
be either, (i) occupied by an agent from subpopulation 1, 
Cw,-,) = 1, (ii) occupied by an agent from subpopulation 2, 

C 2(ij) = l > or N vacant with C T(ij) =0 and C 2(v) = 0 - From our 
simulations we can estimate the average occupancy of agents from 

subpopulation 1 at site (ij) as <Ci(y)> = (l/M) J2^ = l C™ I;) , and 

the average occupancy of agents from subpopulation 2 at site (ij) 

as<C 2( , v) > = (l/M)E^,C- (iV) . 

Results in Fig. 2E-H show <Ci(y)> and <C 2 (y)> associated with 
the simulations in Fig. 2A-D for M = 100. These averaged profiles 
confirm that both subpopulations spread across the domain with 
time and form two moving fronts, one moving in the positive x— 
direction and the other moving in the negative x-direction. The 
averaged density profiles in Fig. 2E-H confirm that both subpop- 
ulations remain well mixed since we have <Ci (y)> ~ <C 2 (y)> at all 
locations and for all time. Results in Fig. 2M-P show <Ci(y))> and 
<C 2 (y)> associated with the simulations in Fig. 2I-L for M= 100. 
These profiles confirm that two moving fronts of cells form with time 
and that one moves in the positive x -direction and the other moving 
in the negative x-direction. The averaged density profiles in 
Fig. 2M-P show that the two subpopulations do not remain well 
mixed since we see that the leading edge of the moving fronts are 
eventually dominated by subpopulation 1 . 

To describe these averaged simulation results using a continuum 
mathematical framework we form two discrete conservation 
statements for <KCi(y)> and <KC 2 (y)>, which describe the the 
change in average occupancy of subpopulation 1 and 2, 
respectively, at site (ij), during the time interval from time t until 
time t + z. The discrete conservation statements are given by 



and 



P 2 (l + P ) 

= ml 4 +PJ <c 2( ,- W) >(i-<c 1M >- 


-<C 2 (ij)}) 


+ P " (1 ~ Px) <C2(, +1 ,/)Xl-<C 1( , v) >- 


-<C 2 (ij)}) 


+ ^p z ) <c2( .._ 1)>(1 _ <cl(i . )> . 


-<C2 M » 


+ 4 <C 2 (y+l)>(l <Cl(, (/ )> 


-<c 2M » 


4 <C 2ft/ )>(l < Q (i+ijj >- 


-<.C 2 ( i+ ij)}) 


4 <C2(y)>(l (ClC-lj-))- 


-<.C 2 (i-ij)}) 


4 <L2(, V )>(1 <Cl(y+l)>- 


"<C 2 (/,/+!)» 


_4xw z ) <C2(i . )>(1 _ <Ci(i ._ i)> . 


-<C 2 (,v-i)» 


p 2 

+ < C 2 (j+ 1 j) >(1 - < Ci (ij) > - < C 2 oj) 


» 


p2 p 

+ -j- < C 2 (i- ij) >(1 - < Ci (ij) > - < C 2 (ij) 


» 


P 2 

+ -j < C 2 (i J+ 1) >(1 - < Ci (ij) > - < C 2 dj) 


» 


pj 

+ < C 2 (ij- 1) >(1 - < Cl (ij) > - < C 2 (ij) 





(4) 



<5<Ci(, v) > = 



< Ci g - 1 j) >( l - < Ci (^ > - < c 2 (ij) » 



+ 



+ 



+ 



Px) <Ci (/+ i,/)>(l - <Ci (,,,)>- <c 2( , v) » 



^(i+p>) 



^(i-p>) 



p 1 



< Cl (ij _ 1) X 1 - < Cl (ij) > - < C 2 (y) » 

< Cl (y + J ) X 1 - ( Cl (ij) ) - < C 2 (ij) » 

< c i (,,,) xi — < Ci (i+ ij) > - < c 2 (,+ ij) » 

< Cl (y-) X 1 - < Cl (i- 1 j) ) - < C 2 (i _ 1 j) )) 
( Cl (y) >( 1 - < Cl (ij + 1) ) - < C 2 (y + 1) )) 
( Cl (ij) >( 1 - < Cl (ij - 1) > - < C 2 (y- 1) >) 



+ ^<C 1( , + li/) >(l-<C 1( y ) >-<C 2( y ) » 
Pi 

+ ^<C 1( ,_ 1 , /) >(l-<C 1( y ) >-<C 2( y ) » 
+ "f ( Cl (ij + 1 ) >( 1 - < Cl (y ) ) - ( C 2 (y) )) 



(3) 



Positive terms on the right of Equation (3) represent events that 
place an agent of subpopulation 1 at site (ij), while the negative 
terms on the right of Equation (3) represent events that remove 
agents of subpopulation 1 from site (ij). A equivalent interpre- 
tation applies to the terms on the right of Equation (4) with respect 
to agents from subpopulation 2. All the terms on the right of 
Equations (3) and (4) involve factors like <C](y)> and 
(1 — <C](y)) — <C 2 (y)» which we interpret the probability that 
site (ij) is occupied by an agent from subpopulation 1, or the 
probability that site (ij) is vacant, respectively. We interpret 
products of these terms as net transition probabilities which means 
that we are making the usual assumption that the occupancy of 
lattice sites are independent [46-48] . As we shall later demonstrate 
(Fig. 2) this assumption appears to be a reasonable for the 
problems we consider here. 

The discrete conservation statements, given by Equation (3) and 
(4) are related to a system of partial differential equations in the 
appropriate limit as A— >0 and T->0 and the averaged data, 
^Ci(y)> and <C 2 (y)> are written in terms of two continuous 
variables C\(x,y,t) and C 2 (x,y,t). To find this relationship we 
expand all terms in Equations (3) and (4) in a truncated Taylor 
series about site (ij), keeping terms up to 0(A 2 ). Dividing the 
resulting expressions by T, we consider the limit as A->0 and T— >0 
simultaneously, with the ratio (A 2 /t) held constant. In the 
continuum limit, the partial differential equations governing 
C\(x,y,t) and C 2 (x,y,t) can be written as 



+ -f <Cl(y-,)Xl - <Cl(y)> - <C 2( y)», 
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3C, d 



ox ox 



dy 



n s/ (| . < 8S 



(5) 



-vi x — [Cia-^-vij, — [Ci(l-S)] + XiCi(l-S), 



Bt dx 



(l-5)^ + C 2 ^ 

ox ox 



+ £>2 



ay By 



- V2 X ^ [C 2 (l - S)] - v 2y ^ [C 2 (l - 5)] + A 2 C 2 (1 - 5), 



where 



To match the averaged discrete simulation data we use the same 
initial condition as in the discrete simulations 



Ci(x,0)= 



and 



C 2 (x,0)= 



0 < x < 12 mm , 

12 < x < 13 mm , 

13 < x < 25 mm , 



0 < x < 12 mm , 

12 < x < 13 mm , 

1 3 < x < 25 mm . 



(8) 



(9) 
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a^o,t->oV 2t 
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a^o.t^oV 2t 



lim 



= lim 

A^0,t->0 



A^0,i^0\ 2t 



lim 

A^O.t^O \ T 



(6) 



and S(x,y,t) = C\(x,y,t) + C2(x,y,t) is the total density [35,39] 

We note that for the special case where the motion is unbiased 
p x = p y = 0, and that both subpopulations are identical with 
D\ = D 2 = D, X\=Xi = X, we can re-write Equation (5) in terms of 
the total population density as 



Tt 



= D 



8 2 S 



8 2 S 



dx 2 ~^ dy 2 



+ XS(l-S), 



(7) 



which is the two-dimensional analogue of the well known Fisher 
Kolmogorov equation [49,50]. This standard reaction diffusion 
model is a particular case of the more general system derived here. 

We note that Equation (5) is written in terms of the two- 
dimensional (x,y) Cartesian plane. If we consider an initial 
condition, C\(x,y,Q) and C 2 (x,7,0), that is independent of the 
vertical coordinate y, and either periodic or reflecting boundary 
conditions are applied on both boundaries parallel to the x 
coordinate, the solution of Equation (5) is independent of y for all 
?>0 and we have C\(x,y,i) = C\(x,i) and C2(x,y,t) = C 2 (x,?) 
[35,39,40]. These initial conditions and boundary conditions are 
relevant when considering an in vitro experiment in a narrow 
channel geometry, such as a scratch assay [2,3,15] or the discrete 
simulations in Fig. 2A-D and Fig. 2I-L. For other types of assays 
where the geometry is genuinely two-dimensional, such as barrier 
assays [45,51-54], we must consider the complete two-dimension- 
al partial differential equations as demonstrated previously in [40] . 

To investigate how the solution of Equation (5) relates to the 
averaged discrete data in Fig. 2E-H and Fig. 2M-P, we solved 
Equation (5) numerically on 0<x<25 mm with reflecting 
boundary conditions for both subpopulations at both boundaries. 



We solve Equations (5) using a non-iterative operator split method 
[55]. To solve the transport terms in Equation (5) we use central 
difference approximation with mesh spacing Sx, and implicit Euler 
stepping with a time step of St. To solve the reaction terms in 
Equation (5) we use a fourth order Runge-Kutta method with time 
step St [56]. 

The numerical solution of Equation (5), with Z>i=Z) 2 and 
X\ = ^2, is superimposed in Fig. 2E-H where we see that C\{x,i) 
and C 2 (x,?) match the averaged discrete data, <Ci(x,?)> and 
■(C 2 (x,?)>, very well at all locations and for all times considered. 
Similarly, the numerical solution of Equation (5) with D\ ^ D2 and 
X\ = X2 are superimposed in Fig. 2M-P where we also see that 
these solutions match the averaged discrete data very well at all 
locations for all times considered. We would like to reiterate here 
that the key result is that when the subpopulations are identical we 
have C\(x,i) = C 2 (x,?) whereas when we consider distinct 
subpopulation with D\ 7^i) 2 , we observe an influence on the 
spatial and temporal organization of the two subpopulations. In 
particular, we see that the more motile subpopulation, C\(x,t), 
dominates the total population at the leading edge whereas the less 
motile subpopulation, C 2 (x,f), is absent from the leading edge. 
Therefore, locally at the invasive front we have S(x,t)x C\(x,i) 
and C 2 (x,?)~0. We note that these kinds of differences, where 
cells at the leading edge of an invasive population are appear to be 
functionally distinct from cells located well behind the leading edge 
of the invasive population have also been observed experimentally 
in various in vivo [19] and in vitro contexts [13]. 

Discussion and Conclusions 

Mathematical and computational modeling has played an 
important role in improving our understanding of collective cell 
spreading in a range of applications [2-5,13,22]. Despite a range 
experimental evidence that suggests otherwise, most mathematical 
models of collective cell behavior make the simplifying assumption 
that the cell motility rate and the cell proliferation rate are 
constants and do not vary amongst the cell population. These 
kinds of simplifying assumptions give rise to mathematical models 
that take the form of reaction diffusion equations with constant cell 
diffusivity [2-5], or discrete random walk models of collective cell 
behavior where isolated individual agents in the system have 
constant rates of motility [13,15,35,39]. 

In this work we have sought to explore the validity of these 
standard assumptions by performing a transwell assay with 3T3 
cells. By stopping the assay after a short period of time we aimed to 
test the hypothesis that those cells amongst the total population 
with high motility rate would move into the lower compartment of 
the transwell faster than those cells amongst the total population 
with a lower motility rate. Indeed our time lapse data suggests that 
those cells that moved into the lower compartment in a short 
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period of time were, on average, approximately three times more 
motile that those cells remaining in the upper compartment of the 
transwell. We illustrated the role that such variability can have by 
applying an existing model of cell migration through a transwell 
[35] which we generalize so that each agent in the simulation can 
have a distinct motility rate and distinct proliferation rate. Taking 
the simplest possible approach where we consider the total 
population to be composed of two subpopulations, we show that 
the mathematical model predicts very different behavior in the 
transwell assay where we account for differences in the motility 
rate between subpopulation 1 and subpopulation 2. 

We also apply our mathematical model to the situations where 
we idealize a total population of cells as two possibly distinct 
subpopulations, to a scratch assay. Our simulations and analysis 
indicate that when we make the standard assumption that both 
subpopulations have identical cell diffusivity {D\=Dt) and 
identical cell proliferation rate (k\ = X2), with the further 
assumption that the initial condition is a well mixed population 
where both subpopulations are present in equal proportions, we 
observe the formation of moving fronts of cells where both 
subpopulations are well mixed throughout. In contrast, if we 
assume that the subpopulations have distinct cell diffusivities 
{D\ > D2) and identical cell proliferation rates (k\ = X2) our 
modeling shows that the moving fronts of cells that forms is very 
different. In this case the two subpopulations do not remain well 
mixed, and instead we observe that the subpopulation with the 
higher diffusivity dominates at the leading edge of the population. 
This idea that the cells at the leading edge of invasive fronts are 
more motile than their counterparts well behind the leading edge 
is consistent with previous experimental observations [13,19]. 

There are several ways that the modeling results can be 
extended. For example, when we considered the scratch assay 
simulations in Fig. 2 we always made the simplifying assumption 
that both subpopulations were initially present in equal propor- 
tions so that we had C\(x,0) = C2(x,0). In the case that we have 
distinct subpopulations with D\>D2, we note that our main 
result, showing that the two subpopulations do not remain well 
mixed after a sufficiently long period of time, also holds when we 
vary the initial condition. For example, our results in Fig. 2 made 
the assumption that the central region of the domain was equally 
composed of both subpopulations, C\(xfi) = C2(x,0) = 0.5 
(Fig. 21). If, instead, we suppose that subpopulation 2 dominates 
initially by setting Ci(x,0) = 0.1 and C2(x,0) = 0.9 in this central 
region, our modeling framework predicts that subpopulation 1, 
with Di>Z>2, will eventually dominate the leading edge of the 
spreading front despite the fact that there is only a small 
proportion of subpopulation 1 present at the beginning of the 
experiment. 

Another simplifying assumption made here is that we supposed 
that the total population of cells in the system could be idealized as 
two subpopulations. This assumption was invoked so that we could 
illustrate our results as simply as possible and we note that our 
discrete modeling framework, outiined in Section 1.2, and the 
associated continuum partial differential equation description, can 
be generalized so that we can consider dividing the total 
population into an arbitrary number of subpopulations. For 
example, if instead of treating the total population in Fig. 2 as two 
subpopulations, we could consider the total population to be 
composed of N>\ subpopulations. Taking the same approach 
leads to the following system of coupled partial differential 
equations 



dC k 8Jx k djy k 



8t 



dx dy 



+ S k , k =1,2,3,..., AT, (10) 



where Jx k is the flux of subpopulation k in the x-direction, Jy k is 
the flux of subpopulation k in the _y-direction and S k is the source 
term for subpopulation k. These terms can be written as 



j Xk = -D k ( 1 - j2 cj e -§- D ^^\J2 c ')+ v ** c * ( 1 - E ( 



Jy k =-D k [l-J2 C,j ^ ~D k C k A \J2 Q ) + v ky C k ( 1 - £ Q 



s k =x k c k i 1 — ^2 Q I , 



(ii) 



where 



D k = lim 



pk A 2 



v kx = lim 



P k m P^ 



A^o,T^o\ 2t /' 



v ky = lim . , 



x k - 



lim 



(12) 



We note that the question of determining the appropriate number 
of subpopulations, N, which accurately reflects the collective 
behavior of a real population of cells is an open question that 
requires further experimental and theoretical investigation. 

Our experimental methods focused on a transwell assay which 
are typically conducted over an interval of time that is much 
shorter than the cell cycle time [35]. As a consequence, our 
experimental methods were aimed at investigating the variability 
of cell motility amongst the population rather than focussing on 
the variability of the cell proliferation rate. To make such an 
assessment, a different kind of experimental system could be 
considered, such as a barrier assay [45,51-53] or a scratch assay 
[2,3,41], which are often conducted for periods of time that are 
longer than the doubling time. We leave such an investigation of 
the role of variations in the proliferation rate of cells for future 
investigation. 

We conclude with a brief discussion about the limitations of our 
mathematical modelling framework, together with a brief discus- 
sion about the suitability of our mathematical modelling frame- 
work for this particular study. One of the limitations of our 
mathematical model is that it neglects to explicitly incorporate any 
details regarding biochemical cell to cell interactions. It is 
important to acknowledge this limitation since biochemical cell 
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to cell interactions are thought to influence collective cell 
behaviour [20,21]. Since the focus of our work is to explore the 
role of variability amongst a population of cells, it is appropriate 
for us to use a mathematical modelling framework that can 
explicitly examine the role of variability rather than a mathemat- 
ical model that incorporates, potentially complicated, biochemical 
cell to cell interactions. Once again, we would like make the point 
that many traditional mathematical models of collective cell 
behaviour treat the motility of cells as a simple constant value 
across a population of cells [2,3, 1 0, 1 1 , 1 3] . Similarly, many 
traditional mathematical models of collective cell behaviour treat 
the proliferation rate of cells as a simple constant value across a 
population of cells [2,3,10,11,13]. In contrast, our experimental 
data showed that measurements of cell diffusivity from a single 
population of cells can lead to a wide range of cell diffusivity 
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